

#pragma once
#include <stdio.h>
#include <stdint.h>
#include <xmmintrin.h>
#include <mmintrin.h>
#include <immintrin.h>
#include <smmintrin.h>
#include <intrin.h>
#include <math.h>

using namespace System;
using namespace RayDen::Library::Core;
using namespace RayDen::Library::Core::Primitives;
using namespace System::Runtime::InteropServices;

namespace CoreCalc {

	class Vec3f;



	public class SSEHelper
	{

	private:
		static void print128_num(__m128i var)
		{
			uint16_t *val = (uint16_t*)&var;//can also use uint32_t instead of 16_t
			printf("Numerical: %i %i %i %i %i %i %i %i \n",
				val[0], val[1], val[2], val[3], val[4], val[5],
				val[6], val[7]);
		}

		static void print128_num(__m128 var)
		{
			float *val = (float*)&var;//can also use uint32_t instead of 16_t
			printf("Numerical: %f %f %f %f \n", var.m128_f32[0], var.m128_f32[1], var.m128_f32[2], var.m128_f32[3]
				);
		}
	public:

		static inline void SSE_Pix(float x, float y, float z, float w, int max, int* dx, int* dy, int* dz)
		{
			//Min(gamma, round(gamma*clamp(i, 0, 1)))
			__m128 g = _mm_set1_ps(max);
			__m128 unit = _mm_set1_ps(1.f);
			__m128 zero = _mm_set1_ps(0.f);

			__m128 c = _mm_mul_ps(_mm_setr_ps(x, y, z, 0), _mm_set1_ps(w));
			//print128_num(c);
			__m128 res = _mm_min_ps(_mm_set1_ps(max - 1), _mm_round_ps(_mm_mul_ps(g, _mm_min_ps(unit, _mm_max_ps(c, zero))), 0));
			__m128i cp = _mm_cvttps_epi32(res);

			*dx = cp.m128i_i32[0];
			*dy = cp.m128i_i32[1];
			*dz = cp.m128i_i32[2];
		}

		static inline void SSE_SphereIntersect(int sphCount, const float* rorg, const float* rdir, const float* sph_d, float* dist)
		{
			float o[4] = { rorg[0], rorg[1], rorg[2], 0.f };
			float d[4] = { rdir[0], rdir[1], rdir[2], 0.f };

			__m128 org = _mm_load_ps(o);
			__m128 dir = _mm_load_ps(d);
			for (int i = 0; i < sphCount*4; i += 4) {
				float ps[4] = { sph_d[i], sph_d[i + 1], sph_d[i + 2], 0.f };
				__m128 p = _mm_load_ps(ps);
				__m128 op = _mm_sub_ps(org, p);
				__m128 b = _mm_dp_ps(dir, op, 0xFF);
				__m128 l2 = _mm_dp_ps(op, op, 0xFF);

				float det = b.m128_f32[0] * b.m128_f32[0] - (l2.m128_f32[0]) + sph_d[i + 3] * sph_d[i + 3];
				if (det < 0.f) {
					dist[i / 4] = 1e10f;
					continue;
				}
				det = sqrtf(det);
				float t1 = b.m128_f32[0] - det;
				float t2 = b.m128_f32[0] + det;
				if (t1 > 1e-1) {
					dist[i / 4] = t1;
				}
				else if (t2 > 1e-1) {
					dist[i / 4] = t2;
				}
				else {
					dist[i / 4] = 1e10;
				}

			}

		}

		// TODO: Add your methods for this class here.
		static void VectMult(float * src1, float * src2, float * dest, int length) {
			const int c_vectorStride = 4;
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_load_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_load_ps(src2 + i); // MOVAPS
				__m128 vresult = _mm_mul_ps(v1, v2); // MULPS
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectMult(float * src1, float src2, float * dest, int length) {
			const int c_vectorStride = 4;
			float f[4] = { src2, src2, src2, src2 };

			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_load_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_loadu_ps(&f[0]); // MOVAPS
				__m128 vresult = _mm_mul_ps(v1, v2); // MULPS
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectAdd(float * src1, float * src2, float * dest, int length) {
			const int c_vectorStride = 4;
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_load_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_load_ps(src2 + i); // MOVAPS
				__m128 vresult = _mm_add_ps(v1, v2); // MULPS
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectAdd(float * src1, float src2, float * dest, int length) {
			const int c_vectorStride = 4;
			float f[4] = { src2, src2, src2, src2 };

			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_load_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_loadu_ps(&f[0]); // MOVAPS
				__m128 vresult = _mm_add_ps(v1, v2); // MULPS
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectMAdd_E(float * src1, float * src2, float * dest, int length) {
			const int c_vectorStride = 4;
			float f[4] = { 0.5f, 0.5f, 0.5f, 0.5f };

			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_loadu_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_loadu_ps(src2 + i); // MOVAPS
				__m128 v4 = _mm_loadu_ps(&f[0]); // MOVAPS
				__m128 v3 = _mm_loadu_ps(dest + i);
				__m128 vresult = _mm_mul_ps(_mm_add_ps(v1, v2), v4);
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectAverage(float * src1, float * src2, float* src3, float * dest, int length) {
			const int c_vectorStride = 4;
			float f[4] = { 1.f / 3.f, 1.f / 3.f, 1.f / 3.f, 1.f / 3.f };

			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_loadu_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_loadu_ps(src2 + i); // MOVAPS
				__m128 v3 = _mm_loadu_ps(src3 + i); // MOVAPS
				__m128 v4 = _mm_loadu_ps(&f[0]);
				__m128 vresult = _mm_mul_ps(_mm_add_ps(_mm_add_ps(v1, v2), v3), v4);
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectAverage(float * src1, float * src2, float* src3, float* src4, float * dest, int length) {
			const int c_vectorStride = 4;
			float f[4] = { 1.f / 4.f, 1.f / 4.f, 1.f / 4.f, 1.f / 4.f };

			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_loadu_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_loadu_ps(src2 + i); // MOVAPS
				__m128 v3 = _mm_loadu_ps(src3 + i); // MOVAPS
				__m128 v4 = _mm_loadu_ps(src4 + i); // MOVAPS
				__m128 v5 = _mm_loadu_ps(&f[0]);
				__m128 vresult = _mm_mul_ps(_mm_add_ps(_mm_add_ps(_mm_add_ps(v1, v2), v3), v4), v5);
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectMAdd(float * src1, float * src2, float * dest, int length) {
			const int c_vectorStride = 4;
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_loadu_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_loadu_ps(src2 + i); // MOVAPS
				__m128 v3 = _mm_loadu_ps(dest + i);
				__m128 vresult = _mm_add_ps(v3, _mm_mul_ps(v1, v2));
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void VectMAdd(float * src1, float  src2, float * dest, int length) {
			const int c_vectorStride = 4;
			float f[4] = { src2, src2, src2, src2 };
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m128 v1 = _mm_loadu_ps(src1 + i); // MOVAPS
				__m128 v2 = _mm_loadu_ps(&f[0]); // MOVAPS
				__m128 v3 = _mm_loadu_ps(dest + i);
				__m128 vresult = _mm_add_ps(v3, _mm_mul_ps(v1, v2));
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static float VectDot(const float * src1, const float * src2, int length) {
			float result = 0.f;
			for (int i = 0; i < length; i += 4) {
				__m128 v1 = _mm_loadu_ps(src1 + i);
				__m128 v2 = _mm_loadu_ps(src2 + i);
				__m128 vresult = _mm_dp_ps(v1, v2, 0xFF);
				result += vresult.m128_f32[0];
			}
			return result;
		};


		static inline float Dot(float srcx, float srcy, float srcz,
			float dstx, float dsty, float dstz) {
			__m128 v1 = _mm_setr_ps(srcx, srcy, srcz, 0.f);
			__m128 v2 = _mm_setr_ps(dstx, dsty, dstz, 0.f);
			__m128 vresult = _mm_dp_ps(v1, v2, 0xFF);
			return  vresult.m128_f32[0];
		};

		static void VectSqrt(float * src1, float * dest, int length) {
			const int c_vectorStride = 4;
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.
				__m128 v1 = _mm_loadu_ps(src1 + i); // MOVAPS
				__m128 vresult = _mm_sqrt_ps(v1); // MULPS
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static void XYZ(const float * spd, const float * xSpec, const float * ySpec, const float * zSpec, float yInt, float * xout, float* yout, float* zout, int length) {
			/*union {
				__m128 v;
				float a[4];
				} ax, ay, az;*/
			const int c_vectorStride = 4;
			const __m128 f_int = _mm_set1_ps(yInt);
			__m128 rsum = _mm_set1_ps(0.0);
			float xx = VectDot(spd, xSpec, length);
			float yy = VectDot(spd, ySpec, length);
			float zz = VectDot(spd, zSpec, length);

			__m128 r_x = _mm_div_ps(_mm_set_ps(xx, yy, zz, 0.f), f_int);


			/*for (int i = 0; i < length; i += c_vectorStride) {
				const __m128 v_spd = _mm_loadu_ps(spd + i);
				const __m128 v2 = _mm_loadu_ps(xSpec + i);
				const __m128 v3 = _mm_loadu_ps(ySpec + i);
				const __m128 v4 = _mm_loadu_ps(zSpec + i);
				const int mask = 0xFF;
				ax.v = _mm_dp_ps(v_spd, v2, mask);
				ay.v = _mm_dp_ps(v_spd, v3, mask);
				az.v = _mm_dp_ps(v_spd, v4, mask);
				X += ax.a[0];
				Y += ay.a[0];
				Z += az.a[0];
				}*/
			//__m128 r = _mm_div_ps(_mm_setr_ps(X, Y, Z, 0.f), f_int);
			//_mm_store_ps(dest, r);

			*xout = r_x.m128_f32[0];
			*yout = r_x.m128_f32[1];
			*zout = r_x.m128_f32[2];

		};


		static void VectClamp(float* src, float *dest, float min, float max, int length)
		{
			const int c_vectorStride = 4;
			float fm[4] = { min, min, min, min };
			float fmx[4] = { max, max, max, max };
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.
				__m128 v1 = _mm_loadu_ps(src + i); // MOVAPS
				__m128 v2 = _mm_load_ps(fm + i); // MOVAPS
				__m128 v3 = _mm_load_ps(fmx + i); // MOVAPS

				__m128 vresult = _mm_max_ss(v1, v3);
				vresult = _mm_min_ss(vresult, v2);
				_mm_storeu_ps(dest + i, vresult); // MOVAPS
			}
		}

		static float Clamp(float val, float minval, float maxval)
		{
			// Branchless SSE clamp.
			// return minss( maxss(val,minval), maxval );

			_mm_store_ss(&val, _mm_min_ss(_mm_max_ss(_mm_set1_ps(val), _mm_set1_ps(minval)), _mm_set1_ps(maxval)));
			return val;
		}
	};


	public class AVXHelper
	{
	public:
		static const int c_vectorStride = 8;

		// TODO: Add your methods for this class here.
		static inline void VectMult(float * src1, float * src2, float * dest, int length) {
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m256 v1 = _mm256_load_ps(src1 + i); // MOVAPS
				__m256 v2 = _mm256_load_ps(src2 + i); // MOVAPS
				__m256 vresult = _mm256_mul_ps(v1, v2); // MULPS
				_mm256_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};


		static void VectMult(float * src1, float src2, float * dest, int length) {
			float f[8] = { src2, src2, src2, src2, src2, src2, src2, src2, };

			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m256 v1 = _mm256_load_ps(src1 + i); // MOVAPS
				__m256 v2 = _mm256_loadu_ps(&f[0]); // MOVAPS
				__m256 vresult = _mm256_mul_ps(v1, v2); // MULPS
				_mm256_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static inline void VectMAdd(float * src1, float * src2, float * dest, int length) {
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m256 v1 = _mm256_loadu_ps(src1 + i); // MOVAPS
				__m256 v2 = _mm256_loadu_ps(src2 + i); // MOVAPS
				__m256 v3 = _mm256_loadu_ps(dest + i);
				__m256 vresult = _mm256_add_ps(v3, _mm256_mul_ps(v1, v2));
				_mm256_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

		static inline void VectMAdd(float * src1, float  src2, float * dest, int length) {
			float f[8] = { src2, src2, src2, src2, src2, src2, src2, src2, };
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.

				__m256 v1 = _mm256_loadu_ps(src1 + i); // MOVAPS
				__m256 v2 = _mm256_loadu_ps(&f[0]); // MOVAPS
				__m256 v3 = _mm256_loadu_ps(dest + i);
				__m256 vresult = _mm256_add_ps(v3, _mm256_mul_ps(v1, v2));
				_mm256_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};
		static float VectDot(const float * src1, const float * src2, int length) {
			float result = 0.f;
			for (int i = 0; i < length; i += 8) {
				__m256 v1 = _mm256_loadu_ps(src1 + i);
				__m256 v2 = _mm256_loadu_ps(src2 + i);
				__m256 vresult = _mm256_dp_ps(v1, v2, 0xFF);
				result += vresult.m256_f32[0];
			}
			return result;
		};
		static inline void VectSqrt(float * src1, float * dest, int length) {
			for (int i = 0; i < length; i += c_vectorStride) {
				// Vector load, multiply, store.
				__m256 v1 = _mm256_loadu_ps(src1 + i); // MOVAPS
				__m256 vresult = _mm256_sqrt_ps(v1); // MULPS
				_mm256_storeu_ps(dest + i, vresult); // MOVAPS
			}
		};

	};



	public ref class SSE
	{
	public:
		static float VectorDot(array<float>^ a, array<float>^ b, int length)
		{
			pin_ptr<float> a_p = &a[0];
			pin_ptr<float> b_p = &b[0];

			return SSEHelper::VectDot(a_p, b_p, length);

		}

		static void ComputePixel(RgbSpectrum% s1, float f, RgbSpectrum% r, array<float>^ gt, int g)
		{
			int dst[3];
			SSEHelper::SSE_Pix(s1.c1, s1.c2, s1.c3, f, g, &dst[0], &dst[1], &dst[2]);
			int c1 = dst[0];
			int c2 = dst[1];
			int c3 = dst[2];

			r.c1 = gt[c1];
			r.c2 = gt[c2];
			r.c3 = gt[c3];
		}

		static void MaddSse(RgbSpectrum% s1, RgbSpectrum % s2, RgbSpectrum% r)
		{

			float src1[4] = { s1.c1, s1.c2, s1.c3, 0.f };
			float src2[4] = { s2.c1, s2.c2, s2.c3, 0.f };
			float dst[4];

			SSEHelper::VectMAdd_E(&src1[0], &src2[0], &dst[0], 4);

			r.c1 = dst[0];
			r.c2 = dst[1];
			r.c3 = dst[2];
		}


		static void AverageSSE(RgbSpectrum% s1, RgbSpectrum % s2, float f, RgbSpectrum% r)
		{
			float src1[4] = { s1.c1, s1.c2, s1.c3, 0.f };
			float src2[4] = { s2.c1, s2.c2, s2.c3, 0.f };
			float dst[4];

			SSEHelper::VectMAdd_E(&src1[0], &src2[0], &dst[0], 4);

			r.c1 = dst[0];
			r.c2 = dst[1];
			r.c3 = dst[2];
		}


		static void AverageSSE(RgbSpectrum% s1, RgbSpectrum % s2, RgbSpectrum % s3, RgbSpectrum% r)
		{
			float src1[4] = { s1.c1, s1.c2, s1.c3, 0.f };
			float src2[4] = { s2.c1, s2.c2, s2.c3, 0.f };
			float src3[4] = { s3.c1, s3.c2, s3.c3, 0.f };

			float dst[4];

			SSEHelper::VectAverage(&src1[0], &src2[0], &src3[0], &dst[0], 4);

			r.c1 = dst[0];
			r.c2 = dst[1];
			r.c3 = dst[2];
		}


		static void AverageSSE(RgbSpectrum% s1, RgbSpectrum % s2, RgbSpectrum % s3, RgbSpectrum % s4, RgbSpectrum% r)
		{
			float src1[4] = { s1.c1, s1.c2, s1.c3, 0.f };
			float src2[4] = { s2.c1, s2.c2, s2.c3, 0.f };
			float src3[4] = { s3.c1, s3.c2, s3.c3, 0.f };
			float src4[4] = { s4.c1, s4.c2, s4.c3, 0.f };

			float dst[4];

			SSEHelper::VectAverage(&src1[0], &src2[0], &src3[0], &src4[0], &dst[0], 4);

			r.c1 = dst[0];
			r.c2 = dst[1];
			r.c3 = dst[2];
		}

		static void ToXYZ(SampledSpectrum%  X, SampledSpectrum% Y, SampledSpectrum% Z, SampledSpectrum% s1, RgbSpectrum% s2)
		{
			float dst[4];
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;

			pin_ptr<float> xp = &X.spectra.FixedElementField;
			float* x_1 = xp;

			pin_ptr<float> yp = &Y.spectra.FixedElementField;
			float* y_1 = yp;

			pin_ptr<float> zp = &Z.spectra.FixedElementField;
			float* z_1 = zp;


			SSEHelper::XYZ(s_1, x_1, y_1, z_1, SampledSpectrum::yint, &dst[0], &dst[1], &dst[2], SampledSpectrum::nSamples);

			s2.c1 = dst[0];
			s2.c2 = dst[1];
			s2.c3 = dst[2];
		}

		static void MulSSE(SampledSpectrum% s1, SampledSpectrum% s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> s2p = &s2.spectra.FixedElementField;
			float* s_2 = s2p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			SSEHelper::VectMult(s_1, s_2, outp, SampledSpectrum::nSamples);
		};

		static void MulSSE(SampledSpectrum% s1, float s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			SSEHelper::VectMult(s_1, s2, outp, SampledSpectrum::nSamples);
		};

		static void AddSSE(SampledSpectrum% s1, SampledSpectrum% s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> s2p = &s2.spectra.FixedElementField;
			float* s_2 = s2p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			SSEHelper::VectAdd(s_1, s_2, outp, SampledSpectrum::nSamples);
		};

		static void AddSSE(SampledSpectrum% s1, float s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			SSEHelper::VectAdd(s_1, s2, outp, SampledSpectrum::nSamples);
		};



		static void VSqrt(Float4% arg, Float4% res)
		{
			pin_ptr<float> pa = &arg.x;
			pin_ptr<float> po = &res.x;
			SSEHelper::VectSqrt(pa, po, 4);
		}

		static void MAddSSE(SampledSpectrum% s1, SampledSpectrum% s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> s2p = &s2.spectra.FixedElementField;
			float* s_2 = s2p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			SSEHelper::VectMAdd(s_1, s_2, outp, SampledSpectrum::nSamples);
		}

		static float Dot(SampledSpectrum% s1, SampledSpectrum% s2)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> s2p = &s2.spectra.FixedElementField;
			float* s_2 = s2p;
			return SSEHelper::VectDot(s_1, s_2, SampledSpectrum::nSamples);
		}

		static float Dot(RgbSpectrum% s1, RgbSpectrum% s2)
		{
			return SSEHelper::Dot(s1.c1, s1.c2, s1.c3,
				s2.c1, s2.c2, s2.c3);
		}

		static float Dot(Vector% s1, Vector% s2)
		{
			return SSEHelper::Dot(s1.x, s1.y, s1.z,
				s2.x, s2.y, s2.z);
		}

		static void MAddSSE(SampledSpectrum% s1, float s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			SSEHelper::VectMAdd(s_1, s2, outp, SampledSpectrum::nSamples);
		}



		static void ClampSSE(SampledSpectrum% s1, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			pin_ptr<float> s2p = &r.spectra.FixedElementField;
			SSEHelper::VectClamp(s1p, s2p, 0.f, 0.99f, SampledSpectrum::nSamples);
		}


		static float Clamp(float val, float minval, float maxval)
		{
			return SSEHelper::Clamp(val, minval, maxval);
		}

		static void SSE_Spheres_Intersect(int count, Point% org, Vector% dir, array<Float4>^ sphs, [Out] array<float>^% dists)
		{
			pin_ptr<float> rorg = &org.x;
			pin_ptr<float> rdir = &dir.x;

			pin_ptr<float> sph = &sphs[0].x;
			float* d = new float[count];
			SSEHelper::SSE_SphereIntersect(count, rorg, rdir, sph, d);
			dists = gcnew array<float>(count);
			for (int i = 0; i < count; i++) {
				dists[i] = d[i];
			}
		}

	};



	public ref class AVX
	{
	public:

		static float Dot(SampledSpectrum% s1, SampledSpectrum% s2)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> s2p = &s2.spectra.FixedElementField;
			float* s_2 = s2p;
			return AVXHelper::VectDot(s_1, s_2, SampledSpectrum::nSamples);
		}
		static void Mul(SampledSpectrum% s1, SampledSpectrum% s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> s2p = &s2.spectra.FixedElementField;
			float* s_2 = s2p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			AVXHelper::VectMult(s_1, s_2, outp, SampledSpectrum::nSamples);
		};

		static void Mul(SampledSpectrum% s1, float s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			AVXHelper::VectMult(s_1, s2, outp, SampledSpectrum::nSamples);
		};

		static void MAdd(SampledSpectrum% s1, SampledSpectrum% s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> s2p = &s2.spectra.FixedElementField;
			float* s_2 = s2p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			AVXHelper::VectMAdd(s_1, s_2, outp, SampledSpectrum::nSamples);
		}

		static void MAdd(SampledSpectrum% s1, float s2, SampledSpectrum% r)
		{
			pin_ptr<float> s1p = &s1.spectra.FixedElementField;
			float* s_1 = s1p;
			pin_ptr<float> resp = &r.spectra.FixedElementField;
			float* outp = resp;
			AVXHelper::VectMAdd(s_1, s2, outp, SampledSpectrum::nSamples);
		}
	};


	public ref class Float4SSE
	{
	public:

	};

}
